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ABSTRACT 

We present the rest-frame V-band luminosity function (LF) of galaxies at 0.4 < z < 4.0, measured from a 
near-infrared selected sample constructed from the NMBS, the FIRES, the FIREWORKS, and the ultra-deep 
NICMOS and WFC3 observations in the HDFN, HUDF, and GOODS-CDFS, all having high-quality optical 
to mid-infrared data. This unique sample combines data from surveys with a large range of depths and areas in 
a self-consistent way, allowing us to (1) minimize the uncertainties due to cosmic variance; and (2) simultane- 
ously constrain the bright and faint ends with unprecedented accuracy over the targeted redshift range, probing 
the LF down to 0.1 L* at z ~ 3.9. We find that (1) the faint end is fairly flat and with a constant slope from 
z = 4, with a = -1.27 ±0.05; (2) the characteristic magnitude has dimmed by 1.3 mag from z ~ 3.7 to z = 0.1; 
(3) the characteristic density has increased by a factor of ~ 8 fromz ~ 3.7 to z = 0.1, with 50% of this increase 
from z ~ 4 to z ~ 1.8; and (4) the luminosity density peaks at z ~ 1 - 1 .5, increasing by a factor of ~ 4 from 
z = 4.0 to z ~ 1 — 1.5, and subsequently decreasing by a factor of ~ 1.5 by z = 0.1. We find no evidence for 
a steepening of the faint-end slope with redshift out to z = 4, in contrast with previous observational claims 
and theoretical predictions. The constant faint-end slope suggests that the efficiency of stellar feedback may 
evolve with redshift. Alternative interpretations are discussed, such as different masses of the halos hosting 
faint galaxies at low and high redshifts and/or environmental effects. 

Subject headings: galaxies: evolution — galaxies: formation — galaxies: fundamental parameters — galaxies: 
high-redshift — galaxies: luminosity function, mass function 



1. INTRODUCTION 

In the current paradigm of structure formation, dark mat- 
ter (DM) halos build up in a hierarchical fashion through the 
dissipationless mechanism of gravitational instability. The as- 
sembly of the stellar content of galaxies is instead governed 
by much more complicated physical processes, often dissipa- 
tive and non-linear, which are generally poorly understood. 
To counter this lack of understanding, simplified analytical 
prescriptions of complex physical processes are employed in 
galaxy formation models. One of the fundamental tools to 
constrain the free parameters of these models is the luminos- 
ity function (LF), as its shape retains the imprint of galaxy 
formation and evolution processes. 

The underlying physical processes that shape the faint end 
of the LF are generally associated with feedback from super- 
novae that is effective in heatin g gas and driving w inds in shal- 
low gravitational potentials (iDekel & Silld I1986L Suppres- 
sion of gas cooling in low-mass halos due to a background of 
photo-ioni zing radiation also contributes in shaping the faint 
end (e.g., Ben son et al.1 12002b . Matching the bright end of 
the LF has proven more challenging. Recent implementation 
of active galactic nucleus (AGN) feedback in semi-analytic 
models has yielded faithful reproductions of the observed 
local rest-frame op ti cal and near-infrar e d global LFs (e.g. 
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Granato etail l2004t iBower et all l2006t ICroton et all 
Menci et all 120061; iKang et al.l 120061; iMonaco et al.1 
Somerville et al.ll2008l) . These m odels provide good m atches 
to the LFs out to z ~ 1.5 (e.g., ISomerville et al.l 120 1 It and 
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references therein), although significant dis agreements with 
obser vations are still present at z > 2 (e.g., iHenriques et al.l 
l20Tl . 

Within the hierarchical structure formation paradigm, one 
naturally expects strong evolution with redshift of the faint- 
end slope with steeper slopes at earlier time, considering 
that the slope of the dark matter mass function is very steep 
(c(dm ~ _ 2), and that the objects that form in these halos con- 
tinue to grow b y continued star formation and mergers with 
each other (e.g., Khochfar & Burkert 2001), hence flattening 
the slope. 

Observationally, evidence of evolution of the faint-end 
slope a has been recently clai med, with a ste e penin g with 
redshift dRyan et al.l 120071) . iKhochfar et aTJ (120071) have 
shown that simulations of galaxy formation and evolution 
are able to reproduce the observed redshift evolution of a. 
Specifically, while supernova feedback is responsible for flat- 
ter a with respect to «dm at any given redshift, its evolution 
with redshift appears to be driven by a steeper dark matter 
mass function at early cosmic times for the range of halo 
masses hosting sub-L * galaxies (e.g., Kh ochfar et al.l 1200/1; 
iKobayashi et afll2007l) . 

Despite the significant improvements in constraining the 
faint-end of t he rest-frame UV LF of dropout galaxies out 
to z ~ 8 (e.g., IOeschetalJl2010bt iBou wens et~alj 120 lU and 
references therein), the faint-end slope of the rest-frame opti- 
cal LFs of galaxies at z > 1 from near-infrared (NIR) selected 
samples, which arguably allow for the assembly of represen- 
tative samples of high-redshift galaxies, is still quite uncer- 
tain. Moreover, the evidence of a steepening of a with red- 
shift comes from several surveys, each with unique selection 
effects, observational biases, and measured at different rest- 
frame wavelengths. Specifically, all measurements at z > 3.5 
come from studies in the rest-frame far-UV, whereas at z < 3 
LFs are typically measured in the rest-frame optical. 
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In this paper, we derive the rest-frame optical (V-band) 
LFs of galaxies over the redshift interval 0.4 < z < 4.0 us- 
ing a NIR-selected composite sample of galaxies constructed 
from the NMBS, the FIRES, the GOODS-CDFS surveys, 
and ultra-deep Hubble Space Telescope (HST) data that self- 
consistently combines the advantages of deep, pencil beam 
surveys with those of shallow, wider surveys. This NIR com- 
posite sample allows us to (1) minimize the uncertainties due 
to cosmic variance, and empirically quantify its contribution 
to the total error budget by exploiting the large number of in- 
dependent field of views; and (2) simultaneously constrain the 
bright and faint ends of the LF with unprecedented accuracy 
and statistics over the entire targeted redshift range. This pa- 
per is structured as follows. In Section|2] we present the com- 
posite NIR-selected sample used to measure the LFs of galax- 
ies at 0.4 < z < 4.0; the methods used to measure the LF (the 
1 /V m ax and the maximum likelihood methods) are presented 
in Section [3] as well as the LFs of galaxies at 0.4 < z < 0.7, 
0.76 <z < 1.1, 1.1 <z < 1.5, 1.5 < z < 2.1, 2.1 <z < 2.7, 
2.7 < z < 3.3, and 3.3 < z < 4.0, and the evolution of the 
Schechter function parameters and luminosity density. Our 
results are summarized in Section [4] We assume Hm = 0.3, 
£1a = 0.7, and Ho = 70 km s _1 Mpc _1 . All magnitudes are in 
the AB system. 

2. THE NIR COMPOSITE SAMPLE 

The data set we use to measure the LF consists of a com- 
posite NIR-selected sample of galaxies built from several 
multi-wavelength surveys, all having high-quality UV to mid- 
infrared photometry: the F aint InfraRed Extragalactic Survey 
(FIRES: iFranx et alJl2003l). the Great Observato ries Origins 
Deep Survey (GOODS: iGiavalisco et alJf2004h. the NEW- 
FIRM Medium-Band Survey (NMBS; Ivan Dokkum et alj 
120091) . the ultra-deep NICMOS observations over the 
HDF-North (HDFN) GOODS field (Dicki nson etaP \l99§L 
Thompson et al. 1999), the ultra-deep WFC3/IR observations 
taken as part of the HUDF09 program (GO 11563; PI: 
Illingworth) over the Hubble Ultra-Deep Field (HUDF), and 
the wide-area WFC3 Early Release Science (ERS) obser- 
vations over the CDF-South GOODS field (GO 11359; PI: 
O'Connell). 

For the FIRES, GOODS-CDFS, and NMBS surveys we 
have adopte d the publi cl y-ava ila ble .g-band s electe d cat- 
alogs from lLabbe et all d2003l) . IWuvts et alj d2008l) . and 
WhitakeretaUMl 3), respectively. We refer to these works 
for a complete description of the observations, reduction pro- 
cedures, and the construction of the ^-selected multiwave- 
length catalogs. Briefly, the FIRES -HDFS catalog (HDFS, 
hereafter) has 833 sources down to Kg 1 = 26.0 over an area of 
2.5' x 2.5'. The FIREWORKS -CDFS (CDFS, hereafter) K s - 
selected catalog comprises 6308 objects down to Kg 1 = 24.6 
over a total surveyed area of 138 arcmin 2 . The NMBS (COS- 
MOS and AEGIS, hereafter) K-selected catalogs comprise 
31,306 and 27,572 objects, respectively, down to Kg 1 = 24.5 
and over a total surveyed area of 0.5 deg 2 . 

For the remaining fields, we have constr ucted //i6o-selecte d 
catalogs following previous studies (e.g., Labbe ^taT]|2003l) . 
The full details of the reduction, source detection, and genera- 
tion of the photometric catalogs will be described in Stefanon 
et al. (2012, in prep.). Briefly, simil arly to the other adopted 
catalogs, we have used SExtractor (Be rtin & Arnoutslll996l) 
in dual-image mode, with fluxes measured in the registered 
and PSF-matched images in order to derive accurate colors 
and limit any bandpass-dependent effects. We have used an 



alternative source fitting algorithm especially suited for heav- 
ily confused images for which a higher resolution prior (in 
this case, the //)6o-band image) is available to extract the pho- 
tometry from the IRAC images. This method is describe d 
in more detail in the appendix of iMarchesini et ail d2009l) . 
The UDF WFC3 (UDF, hereafter), GOODS-CDFS WFC3 
(ERS, hereafter), and the HDFN NICMOS (HDFN, hereafter) 
Hi6o-selected catalogs comprise 2100, 9470, and 1774 objects 
down to H\(,o = 28.5, Hieo = 27.3, and H\so = 28.0, over an area 
of 5.3, 46.2, and 6.7 arcmin 2 , respectively. 

Table Q] provides a summary of the catalogs used, includ- 
ing the filters that each field has been imaged with, the ob- 
served K- and H-band limiting magnitudes adopted to con- 
struct the composite sample {Kg^ and Hg^; see § 12.11 ), the 
effective area of each pointing (A e ff), the number of galax- 
ies over the redshift interval 0.4 < z < 4.0 used to measure 
the rest-frame V-band LFs (-/Vo.4< z <4.oX an d the correspond- 
ing fraction of galaxies with spectroscopic redshift (/ zspe c)- 
For the HST Z/i6o-selected catalogs, A e ff represents the area 
in the fields where the NIR observations are at least 10% of 
their maximum depths in the field, and with full coverage in 
the ACS bands. For the HDFS and FIREWORKS catalogs, 
A e ff represents the area in the fields over which the optical and 
NIR observations are at least 20% of their maximum depths in 
the field, follow i ng s uggested criteria by Labbe et alj d2003l) 
and lWuvts etaLldioOSfe . 5 Finally, for the NMBS catalogs, A eB 
represents the area in the fields over which the optical and 
NIR observations are at least 30% of their max imum depths in 
the fie ld, again following suggested criteria by Whi tak er et alj 
d20llh . The regions around bright stars were also excluded in 
the calculation of the LF, as faint galaxies in these regions are 
either obscured by the foreground star or have systematically 
incorrect magnitudes due to scattered light. 

2.1. Completeness 

The ffi60- and ^T-band limiting magnitudes listed in Ta- 
ble Q] were chosen following a conservative approach that en- 
sures completeness better than 90% in all catalogs down to 
the adopted limiting magnitudes. The completenesses for the 
FIRE S and F IREWORK S surveys have bee n characterized in 
lLabbe et alj d2003l) and IWuvts et alj d2008l) . For these cata- 
logs, we have used the same /T-band limitin g magnitudes pre- 
viously adopted in Marc hesini et alj d2009h . 

The c ompletenesses of the N MBS catalogs have been de- 
rived in Whita ker et alj d201 lb . When masking sources, the 
point-source completeness is better than 90% at Kg 1 < 23.2. 
When objects are not masked, objects can fall on or close to 
other sources, and these blended objects are not properly han- 
dled by SExtractor. In this case, the completeness in COS- 
MOS and AEGIS is better than 85 and 87% at Kg* < 22.7 
and 22.8, respectively. However, in our work we have adopted 
the "de-blended" catalogs, for which these comp leteness lim- 
its are quite conservative (Whitak er et al.l [20 1 ll) . To derive 
robust limiting magnitudes for the NMBS de-blended cata- 
logs used in our work, we have used the raw number counts 
of galaxies, which do not show any signs of incompleteness 
down to Kg 1 « 22.8, implying completeness better than 95% 
in the de-blended NMBS catalogs at Kg* < 22.7 and 22.8 
for COSMOS and AEGIS, respectively. Therefore, we have 
adopted Kg* = 22.7 and Kg* = 22.8 for COSMOS and AEGIS, 
respectively, as the limiting magnitudes to construct the NIR- 

5 The area of FIREWORKS listed in TablefT]excludes the area covered by 
the ERS and the UDF catalogs. 
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FIG. 1 . — Spectroscopic vs. photometric redshifts for the NIR-selected cat- 
alogs, shown on a "pseudo-log" scale. The number and fraction of sources 
with spectroscopic redshifts used in the shown comparison are also specified, 
as well as the median in Az/(1 + z S pec), ffNMAD. and the fraction of catas- 
trophic outliers. The comparison between Zspcc and z p hoi is extremely good, 
both at low and high redshifts. 



selected composite sample. 

The completenesses for the HST i/ieo-selected catalogs 
(UDF, HDFN, and ERS) will be presented in detail in Ste- 
fanon et al. (2012; in prep.). Briefly, the comp leteness in each 
catalog was derived following the method in Whitak er et al.1 
d201 ll) The (point-source) completeness was estimated as 
a function of magnitude by attempting to recover simu- 
lated sources within the Z/ieo-band noise-equalized detec- 
tion image. The derived completeness is better than 90% at 
H m <28.2, 26.7, and 26.6 for the UDF, ERS, and HDFN 
catalogs, respectively. We have therefore adopted the conser- 
vative limiting magnitudes listed in Table Q] to construct the 
NIR-selected composite sample. 

2.2. Photometric Redshifts 

The overall fraction of sources with spectroscopic redshifts 
in the i/i6o- and ^-selected catalogs is 5.4%. Consequently, 
we must rely primarily on photometric redshift estimates. 
Photometric redshifts z p h t for all galaxi es were derived usin g 
the EAzY photometric redshift code {Brammer et al. 2008), 
adopting z_peak as the photometric redshift, which finds dis- 
crete peaks in the redshift probability function and returns the 
peak with the largest integrated probability. The temp late set 
used in this work is described in Whita ker et al.1 (1201 ll) . com- 
posed of EAzY default templates and an additional template 
for an old, red galaxy. We quantify the accuracy of the photo- 
metric redshifts z p h t by comparing them to the available spec- 
troscopic redshifts z spec , as shown in Figure Q] 

The spectroscopic redshifts for the NMBS COSMOS and 
AEGIS catalogs were taken from the zCOSMOS survey 
(iLillv et al.l 120071) and the Deep Extragalactic Evolutionary 
Probe 2 dDavis et al.ll2003l) . respectively. Spectroscopic red- 
shifts of 1 17 Lyman break g alaxies at z ~ 3 with in the AEGIS 
field were also added from lSteidel etal.1 (120031) . For FIRES 
and FI REWORKS we used the s pectroscopic redshi fts com- 
piled in lRudnick et al.l(l2003l) and lWuvts etaTI d2008l) , respec- 
tively. For the UDF and ERS //[6o-selected catalogs, spec- 
troscopic redshifts were taken from the "CDFS Master Cat- 
alog" of spectroscopic redshifts (v2.0; I. Balestra; Decem- 




FlG. 2. — Left: SEDs of galaxies at 2.1 < z < 2.7 taken from the seven 
NIR-selected catalogs used to construct the composite sample: ERS, HDFN, 
FIREWORKS CDFS, NMBS COSMOS, NMBS AEGIS, FIRES HDFS, and 
UDF (from top to bottom). The galaxies were arbitrarily selected to have 
photometric redshifts falling in the targeted redshift interval, with reason- 
ably well-behaved redshift probability functions. They are not representative 
of the galaxy population at 2.1 < z < 2.7. Circles are the observed fluxes 
in units of F\, with 1 a errors. The solid gray curves represent the best-fit 
EAzY templates; each SED was normalized and offset with respect to the oth- 
ers.Right: EAzY redshift probability functions (gray regions). The red line 
is the adopted redshift from EAzY (z pe ak); the blue line is the spectroscopic 
redshift, when available. The yellow regions are the 1 <r allowed values for 
Zpeak- Also listed are the K- or //-band total magnitudes, z pe ak> and Zspcc. 



ber 2009). 6 The spectroscopic r edshifts for the HDF N Z/ieo- 
selected catalog were taken from lBarger et al.l d2008l) . 

The comparison between spectroscopic and photometric 
redshifts results in excellent z p hot estimates. The median in 
Az/(1 +z sp ec), with Az = Zphot-Zspec, is -0.000, -0.002, -0.007, 
0.018, -0.005, -0.009, and -0.015 in COSMOS, AEGIS, 
CDFS, HDFS, HDFN, ERS, and UDF, respectively, with the 
normalized median absolute deviation (7nmad 7 =0.008, 0.015, 
0.029, 0.040, 0.029, 0.037, 0.039. The fraction of catastrophic 
outliers (defined as galaxies with Az/(l+z spec ) > 5onmad) * s 
4.8%, 3.2%, 2.9%, 1.9%, 2.4%, 3.6%, and 4.6% in COSMOS, 
AEGIS, CDFS, HDFS, HDFN, ERS, and UDF, respectively. 
Figure [2] provides examples of spectral energy distributions 
(SEDs) of galaxies at 2.2 < z < 2.7 in the seven NIR-selected 
catalogs, along with their corresponding EAzY redshift prob- 
ability functions. 

The accuracy and reliability of the photometric red- 
shifts estim ated for the NMBS su rvey have been discussed 
in detail in Whita ker et alj d201 ll) . where excellent agree- 
ment with the available spectroscopic redshifts is found. 
The accuracy of the photometric redshifts estimated in the 
FIRES a nd FIREWORK S catal o gs have been discus sed in 
detail in iBrammer et al] d2008l) . IWuvts et alj d2008l) . and 
iMarchesini et alj d2009l) . A detailed discussion of the qual- 

6 |http : / /www, eso ■ org/ s ci/ activities /pro jects /goods /Master Spec 
7 The normalized median absolute deviation unmad, defined as 1.48 X 
median[\(Az — median(Az))/(i +z S pec)|], is equal to the standard deviation 
for a Gaussian distribution, and it is les s sensitive to outliers than the usual 
definition of the standard deviation (e.g. Brammer et al. 2008) 
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Table 1. Summary of catalogs used to construct the NIR-selected com- 
posite SAMPLE 



Field/Survey 


Selection Band 


Filter Coverage 


K/H! ot 

' lim 
(mag) 


A e ff 
(arcmin 2 ) 


No.4<z<4.0 


/zspec 
[%] 


Ref. 


UDF 


WFC3 Hi«i 


ACS B/S^Vfin^lTT^/ein 

NICMOS J] 10H160 
WFC3 Y 105 J 125 H, 6 o 
ISAAC K s , IRAC 


27.8 


5.1 


547 


7.9 


(1,2,3,4) 


ERS 


WFC3 H 160 


ACS B43 5 V 6 06l775 z 850 

WFC3 UV 2 25,275.336Yl05 T 125Hl60 

ISAAC HK S , IRAC 


26.3 


42.5 


2437 


10.1 


(3,4,5,6,7) 


HDFN 


NICMOS H 160 


AC ,N H -! 3^ V &f\&]T7< 7,QK(\ 

WFPC2 U TO oB 4 5oV 6 o6l8l4 
NICMOS J, i H 160 , IRAC 


26.1 


6.7 


435 


15.6 


(4,5,6,8,9,10) 


HDFS/FIRES 


ISAAC K s 


WFPC2U300B450V6O6I814 
ISAAC J S HK S , IRAC 


25.6 


4.5 


204 


8.3 


(11,12,13) 


CDFS/FIREWORKS 


ISAAC Ks 


ACS B 43 5V 6 06i775 z 850 
WFI U 38 BVRI, IRAC 
ISAAC J S HK S , IRAC CH1:4 


23.2 


77.2 


689 


44.6 


(5,14) 


AEGIS/NMBS 


NEWFIRM K 


GALEX FUV 15 ooNUV 2 300 
CFHTLS ugriz 
NEWFIRM J] J?J}HiHiK 
WIRCam JHK S , IRAC 


22.8 


743 


7402 


15.9 


(15,16) 


COSMOS/NMBS 


NEWFIRM K 


GALEX FUV 15 ooNUVi 30 o 


22.7 


741 


7689 


7.3 


(15,16) 



Subaru BjVjr'i'z', CFHTLS ugriz 
Subaru 12 medium-bands 
NEWFIRM J, J 2 J 3 H,HoK 
WIRCam JHK S , IRAC 



NOTE. — ^V^ihi ' s tne observed K- and //[go-band limiting magnitudes adopted to construct the NIR-selected composite sample; A c ff is the effective area; 
Afo.4<z<4.0 is the number of galaxies over the redshift interval 0.4 < z < 4.0 used to measu re the rest-fr ame V-band LFs, and f z s D e c is t he corresponding frac tion 
of galaxies with z S pec ■ The referenc es of the used catalogs a nd data sets are (DIThompson et alj <200 5). (2) Beckwith et al. 12006), (3) Retzla ffetal] BOlOl). (4) 
Stefan on et al. (2012; in prep ), (5)lGiavalisco et a l. 12004), (6) Dickinson et al. (2003), (7) Damen et al. (201 1), (8) Will iams et alj fl 996). (9) Dickin son et all 

(ll) IFranx et al.N2003l) . (12) ILabbe et alj 120031) . (13) IWuvts etaiT)2007l) , (14) IWuvts et al.N2008l) . (15) lvan Dokkum et all 



1999), (10) 



2009), (16) 



Thompson 



Whitaker 



in prep. ), p) |y 
iret^njl999l) , ( 
et al] 1201 II). 



ity of the photometric redshifts in the Z/160-selected catalogs 
in UDF, ERS, and HDFN will be presented in Stefanon et al. 
(2012; in prep.). 

It should be stressed that the assessment of the errors of 
the photometric redshifts from the z p h t versus z spec compari- 
son works only if the spectroscopic sample is a random sub- 
set of the photometric sample. This is not true in general, 
as the currently available spectroscopic samples are heav- 
ily weighted to ward bright and/or bl ue, low-redshift galax- 
ies. However, [Brammer et al. (2008) showed explicitly that 
the errors in z p hot versus z spec plots are generally consistent 
with the errors estimated from both Monte Carlo simulations 
and from the formal errors in the fit, i.e., the confidence inter- 
vals output by EAzY correctly describe the deviations from 
the spectroscopic redshifts. We are therefore confident that 
the photometric redshift uncertainties returned by EAzY are 
reliable. 

2.3. Rest-frame Luminosities 

Rest-frame luminosities are derived by integrating the red- 
shifted rest-frame filte r bandpass from t he bes t-fit EAzY tem- 
plate, as described in Bram mer et al.l (1201 ll) . This method 
produces very similar rest-frame luminosities as other meth- 
ods interpolating between the observed bands that bracket th e 
rest-frame band at a given redshift (e.g., Rud nick et al.ll2003l) . 
Specifically, down to the H160- and /T-band limiting magni- 
tudes of the NIR-selected composite sample, the difference 
between the rest-frame V-band magnitudes derived with the 
two methods are very small, ~ 0.02 mag, and no dependence 
on the //160- or /T-band magnitude. We note that most galaxies 
in the NIR-selected composite sample have robust photome- 



try in the IRAC bands , which allows us to properly model 
the rest-frame optical SEDs of galaxies and to robustly derive 
their rest-frame V-band magnitudes out to the highest targeted 
redshift. We compu ted rest-frame lumin osities in the V pass- 
band defined by Mafz Apellaniz ( 2006). In all cases where 
a spectroscopic redshift is available, we computed the rest- 
frame luminosity fixed at z spe c- 

2.4. The Composite Sample 

We constructed a composite NIR-selected sample of galax- 
ies over the redshift range 0.4 < z < 4.0 to be used in the 
derivation of the rest-frame V-band LFs of galaxies in § [3] 
The final composite sample includes 19403 NIR-selected 
galaxies with Hf£ Q < 27.8 or K tot < 25.6, and 0.4 < z < 4.0 
over an effective area of 1620 arcmin 2 . Of these sources, 
12.5% have spectroscopic redshifts. This NIR-selected sam- 
ple is unique in that it combines data from surveys with a large 
range of depths and areas in a self-consistent way. The wide, 
fairly deep, NMBS survey, with its exquisite photometric red- 
shifts, provides optimal sampling of the bright end of the LF 
down to the characteristic magnitude at z ~ 3 and a factor of 
~ 20 fainter than the characteristic luminosity, L*, at z ~ 0.6. 
The ultra-deep i/160-selected catalogs allowed us to probe the 
faint end down to 0.1 L* at z ~ 3.9, and to 0.003 L* at z ~ 1 
(see also Figure[3]l. 

3. THE LUMINOSITY FUNCTION AND DENSITY 

8 For the NIR-selected composite sample, 4.1% of the sources at 0.4 < 
z < 4.0 have IRAC signal-to-noise ratio (S/N) <3 in the 3.6 \xm band. This 
fraction is only slightly larger (5.5%) in the highest targeted redshift interval 
3.3 < z < 4.0. 
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3.1. Methodology and Uncertainties 

To estimate the rest-frame V-band LF in the case of a 
composite sample, we have us ed the extended version of the 
1/Vmax estimator as defined in lAvni & Bahcalll (Tl980). The 
Poisson error in each magnit ude bin was computed adopting 
the recipe of iGehrel i d 1986b . The 1/V max estimator has the 
advantages of simplicity and no a priori assumption of a func- 
tional form for the LF; it also yields a fully normalized solu- 
tion. However, it can be affected by the presence of clustering 
in the sample. Uncertainties due to field-to-field variations in 
the determination of the LF with the 1 /V max method were esti- 
mated following the procedure presented in Marchesi ni et al.1 
d2007l) . and added in quadrature to the Poisson errors. The 
contribution from cosmic variance to the error budget is typ- 
ically 0.06 dex at the characteristic magnitude in all targeted 
redshift intervals, except at 3.3 < z < 4.0, where it increases 
to 0. 1 dex. 

We also measured the LF u sing the STY method 
(ISandage. Tammann & Yahill 1 19791) . which is a paramet- 
ric maximum likelihood estimator. The STY method is 
unbiased with respect to density inhomogeneities (e.g., 
lEfstathiou. Ellis & Peterson! 1 19881). it has well-de fined 
asymptotic error properties (e.g jKendail & Stuartlll96il) . and 
it does not require binning of the data. For the ST Y method, 
we ha ve assumed that the LF is described by a ISchechterl 
(119761) function, 

$(M) = (0.4 In 10)$* [io°- 4 < M *- M X I+ «>] 

■exp[-10 a4(M *- M) ], (1) 

where a is the faint-end slope parameter, M* is the charac- 
teristic absolute magnitude, and <£>* is the normalization. The 
best-fit solution is obtained by maximizing the likelihood A 
with respect to the parameters a and M* . The value of $* is 
then obtained by imposing a normalization on the best-fit LF 
such that the total number of observed galaxies in the sam- 
ple is reproduced. We have added in quadrature the contribu- 
tion from cosmic variance previously derived (i.e., 0.06 dex 
at z < 3.3 and 0. 1 dex at 3.3 < z < 4.0) to the error budget of 

The uncertainties on the LF due to photometric red- 
shift errors have been estim ated following the recipe in 
iMarchesini et al.1 (12009, 2010) . Briefly, for each galaxy in the 
NIR-selected sample, a set of 200 mock SEDs was created by 
perturbing each flux point according to its formal error bar. 
Second, photometric redshifts and rest-frame V-band lumi- 
nosities were estimated as described in § [2] Finally, the LFs 
were derived with the 1 /V max and the STY method for each of 
the 200 Monte Carlo realizations of the NIR-selected sample. 
This approach naturally addresses the fact that fainter sources 
tend to be characterized by less accurate z p h t estimates due to 
the larger errors in their photometry, as well as sources charac- 
terized by power-law SEDs and consequently by very poorly 
constrained z p h t estimates and very broad z p h t distributions 
derived from the Monte Carlo realizations. 9 

9 We note that catastrophic outliers in the photometric redshift distribu- 
tion could potent i ally cause systematic err ors in the LF measurements (e.g., 
Chen et al. 2003). Marchesini et al. 12007) performed extensive simulations 
to assess the impact of catastrophic outliers on the derived Schechter func- 
tion parameters, showing that catastrophic outliers will result in measured 
LFs that are typically steeper than the true LFs, and with brighter charac- 
teristic luminosities. For example, at z > 2, they found Aa = —0.01 and 
AM* = —0.05 mag when the effect of 5% catastrophic outliers is simulated, 



The rest-frame V-band LFs of galaxies were derived in the 
following redshift intervals: 0.4 < z < 0.7, 0.76 < z < 1.1, 
1.1 < z < 1.5, 1.5 <z< 2.1, 2.1 <z<2.7,2.7<z<3.3, and 
3.3 < z < 4.0. For all but the lowest redshift interval, the full 
NIR-selected composite sample was used. In the redshift in- 
terval 0.4 < z < 0.7, only the NMBS ^-selected samples were 
used, as the volumes probed by the other surveys become too 
small to be meaningful. We have excluded the redshift inter- 
val 0.7 < z < 0.76, as known clusters of galaxies are present 
in the COSMOS and CDFS fields (e.g., IGuzzo et ail 120071: 
iTrevese et al.ll2007l: rSalimbeni et al. 2009J), contaminating the 
field galaxy LFs presented in this paper. 

Because of the coupling between the two parameters a and 
M*, the luminosity density (obtained by integrating the LF 
over all magnitudes) is a robust way to characterize the evo- 
lution of the LF with redshift. The luminosity density was 
calculated using 

p L = L„<S>(L v )dL„ = T(2 + aWL*, (2) 
Jo 

which assumes that the Schechter parameterization of the ob- 
served LF is a good approximation and valid also at luminosi- 
ties fainter than probed by our composite sample. The uncer- 
tainties on the lumin osity density wer e calcu lated using the 
approach adopted in Marche sini et al.1 (120091) . Specifically, 
the 1 a errors on the total luminosity density have been es- 
timated by deriving the distribution of all of the values of Pl 
allowed within the 1 a solutions of the Schechter function 
parameters from the maximum likelihood analysis. The con- 
tributions to the total error budget from photometric redshift 
random uncertainties (derived with the Monte Carlo simula- 
tions) and from cosmic variance (i.e., 0.06 dex at z < 3.3 and 
0.1 dex at 3.3 < z < 4.0) were also added in quadrature. 

3.2. Results 

The rest-frame V-band LFs of galaxies at 0.4 < z < 0.7, 
0.76 < z < 1.1, 1.1 < z < 1.5, 1.5 < z < 2.1, 2.1 < z < 
2.7, 2.7 < z < 3.3, and 3.3 < z < 4.0 are shown in Fig- 
ure [3] Also plotted are the LFs in the rest-frame r band 
(similar to our rest-frame V band) at z ~ 0.1 from the 
Sloan Digital Sky Survey Six th Data Release (SDSS DR6; 
iMontero-Dorta & Pradall2009l) and the Galaxy and Mass As- 
sembly (GAMA) survey (lLoveday et al.ll2l)lll) . to show the 
evolution from z = 0.55 down to z ~ 0.1. Not e that the local 
LFs fr om IMontero-Dorta & Pr ada (2009|) and lLoveday et alj 
(1201 ll) are essentially identical. Figure |4] shows the evolution 
with redshift of the best-fit values and the 1 a and 2a confi- 
dence contour levels of the two Schechter function parameters 
a and M*, derived with the maximum likelihood analysis. 

The large surveyed area provided by the NMBS catalogs 
allows for the determination of the bright end of the rest-frame 
optical LFs of galaxies at z > 1 with unprecedented accuracy, 
while the depth of the //)6o-selected catalogs in UDF, HDFN, 
and ERS allowed us to constrain the faint-end slope out to 
z ~ 3.5, representing a significant improvement with respect 
to previous works. Specifically, the NMBS samples the LF 
down to the characteristic magnitude at z ~ 3, and down to ~ 

much smaller than the random uncertainties on the Schechter function param- 
eters. Whereas it is likely that these systematics effects get larger at lower 
redshift, especially for the characteristic luminosity, by scaling the results 
from the simulations in Marchesini et al. 1 2007), we expect these systematic 
effects to be in general smaller than, or at most comparable to, the derived 
random uncertainties on the Schechter function parameters at z < 1 ■ 
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FIG. 3. — Rest-frame V-band LFs of galaxies at 0.4 < z < 4.0. Black filled circles are the 1 /V m ax method estimates with total 1-cr error bars. The black solid 
curves are the LFs estimated with the maximum-likelihood method. The shaded gray regions represent the 1-, 2-, and 3-cr uncertainties of the LF estimated from 
the maximum-likelihood method. The blue and red dot-dashed curves show the LFs measured in the adjacent lower and higher redshift interval, to highlight the 
evolution with red shift of the LF. The dot-d ashed cyan curve represents the local (z ~ 0.1) LFs in the rest-frame 'r band (similar to our rest-frame V band) 
from the SDSS I Montero-Dorta & Prada 2009) and the GAMA survey lLovedav et al. 2011). The vertical arrows indicate the approximate 90% completenesses 
in rest-frame V-band magnitude of the catalogs used in our work. 



0.05 L* at z ~ 0.6, while the ultra-deep Hi6o-selected catalogs 
allowed us to probe the faint end down to ~ 0.1 L* at z ~ 3.9, 
and to ~ 0.003 L* at z ~ 1. The approximate completenesses 
in rest-frame V-band magnitude of the catalogs used in our 
work are indicated in Figure|3]with arrows. Most importantly, 
the constructed NIR composite sample allowed us to measure, 
for the first time, the rest-frame optical LF of galaxies from 
z = 4.0 using a single, self-consistent NIR composite sample 
that fully probes both the bright and the faint ends of the LF 
over the entire targeted redshift range 0.4 <z< 4.0. 

Table |2]lists the derived Schechter function parameters and 
the estimated luminosity densities in the targeted redshift in- 
tervals. Figure [5] shows the evolution with redshift of the 
Schechter function parameters a, My, and $*, as well as the 
evolution with redshift of the luminosity density. Figure [5] 
also shows the rest-frame V-band Schechter function param- 
eters and luminosity densities derived i n other works from 
the literatu re, i.e., the Century Survey dBrown et al.l 1200 fl) . 
the SDSS (B an ton et al.l l2003). the VIMOS-VLT Deep Sur- 



vey (VVDS: lllbert et alJl2005T). th e MUSYC-FIREWORKS- 
FIRES surveys (Marchesini et al. 2007), the Cosmic Evo- 
lution Survey (COSMOS ; iLiu et al.l l2008h. the SDSS-DR6 
dMontero-Dorta & Pr ada 2009), and the most recent GAMA 
survey (lLovedav et alJ[20Tlh . 



The faint-end slope a appears to evolve little with redshift. 
The weighted average of our measurements and the z ~ 0.1 
value from the GAMA survey (identical to the SDSS DR6) 
is a = -1.27 ±0.05, shown in Figure[5]by the solid blue line 
and the shad ed gray reg i on. Adopting the simple linear fit 
presented in Rvan et al. (20071) (a(z) = a + bz), we find a = 
-1 .26 ±0.02 and b = -0.02±0.03 (plotted as a black solid line 
in Figure [5]l us i ng the measurements from our work and from 
lLovedav et al.1 (1201 ll) in the fitting analysis. This is consistent 
with a constant faint-end slope in the rest-frame V band over 
the redshift range < z < 4. 10 

We have also modeled the evolution with redshift of the 
characteristic magnitude My and th e normalization <i>* us- 
ing the parameterization presented in Stefanon & Mar chesinll 
(2011). Specifically, the following parameterization was 
adopted to model the observed evolution of M v : 

Mv(z) = (J, [(1 +z)/(l + z*)] "exp [-(1 + Z )/(1 + Z*)] , (3) 

with /i, z*, and 77 free parameters to be determined. By 
performing a least-square fit to our measurements and from 

10 If only the points from our work are used in the fitting analysis, we 
obtain a = —1 .30 ± 0.06 and b = 0.00 ± 0.05, with a resulting weighted mean 
a = —1.29 ±0.05, quantitatively consistent with the analysis including the 
z ~ 0.1 value. 
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Table 2. Best-fit Schechter Function Parameters and Luminosity Densi- 
ties 



Redshift 


My 


a 




log PL 




(mag) 




(10 -4 MpcT 3 mag 


_1 ) (erg sT l Mpc' 3 Hz" 1 ) 


0.4 < z < 0.7 
0.76 < z< 1.1 
1.1 <z.< 1.5 
1.5 <z< 2.1 


-21 76+°' 13 

-Q-H 
97 1^+0.08 

-8-88 

-22 ?0 +8 - i " 

.97 4 7 +0.07 

-Q-O! 


-1.25±0.07 
-1.33±0.04 
-1.28±0.05 
-1.23±0.06 


25.61+ 5 ' 64 
26 lfi**^ 
22.25 + |™ 
18.05|J 


-QQ7. 
96 qq+0-07 
zo - yy -0.08 


2.1 <z< 2.7 
2.7 <z.< 3.3 
3.3 < z< 4.0 


.97 cn+O.II 
J7 -0 12 
-22.83+0-'* 
-22 76+o :40 


-1.37±0.11 
-1.35±0.20 
-1 69+° 58 


9 82+ 2 ' 28 
6 99+2.12 

4 04+4 60 


26.83« 

26 72+0- 16 
zo - /z -0.18 

26.8+};| 



NOTE. — The quoted errors correspond to the 1 a errors estimated from the maximum likelihood analysis, including the contribution from cosmic variance 
and from photometric redshift uncertainties. 
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FIG. 4. — (a-My) parameter space derived from the maximum likeli- 
hood analysis. Filled circles are the best-fit values of a and My derived in 
our work, while the curves represent their lcr and 2cr contour levels. The 
black fi lled square represents the r edshift z ~ 0.1 value from the SDSS- 
DR7 I Montero-Dorta & Prada 2009), whereas the empty gray circle repre- 
sents the redshift z ~ 0. 1 value from the GAMA survey ILoved av et alJ201U) . 
The dotted line represents the weighted average of the plotted measurements 
(a = 1.27 ±0.05). 



ILovedav et alj d201 ll) . we obtain the following values: fx = 
-28.8 ±0.4 mag, z* = 191 ± 338, and rj = (55.5 ± 17.2) x 10" 3 . 
The resulting curve is plotted as a solid dark gray curve in Fig- 
ure [5] The characteristic magnitude is seen to brighten with 
redshift, from M v = -21 .5 mag at z ~ 0.1 to = -22.8 mag 
at z~ 3.7. 

The following parameterization was adopted to model the 
observed evolution of $*: 

<^(z) = 9exp[ 1 /(l+zf], (4) 

where 9, 7, and f3 are the free parameters. The best-fit values 
are 6 = (3.2 ± 0.4) x 10" 3 mag" 1 Mpc" 3 , 7 = (-4.3 ± 5.0) x 
10" 2 , and /3 = -2.65 ± 0.88. The best-fit of Equation|4]is plot- 
ted as a solid dark gray curve in Figure [5] The normalization 
is seen to increase with cosmic time, by a factor of ~ 6^ 2 
from z ~ 3.7 to z ~ 0.6, and by a factor of ~ 8^ 3 to z ~ 0.1, 
with 50% of such increase in the 2 Gyr from z = 4 to z = 1.8, 



and the remaining increase in the following 9 Gyr from z = 1 -8 
toz = 0.1. 

Using the expressions of Eq. [4] and Eq. [3]in Eq. [2] we can 
derive a functional representation of the luminosity density. 
The dashed dark gray curve in the bottom panel of Figure [5] 
represents the luminosity density obtained with this method 
and adopting the best-fit values of the parameters previously 
recovered. The agreement with the points is good over the 
entire redshift range. We stress that no fitting was performed 
using the luminosity density estimates. The luminosity den- 
sity appears to peak at z ~ 1-1.5, increasing by a factor of 
~ 4 from z ~ 3.7 to z ~ 1 - 1 .5, and subsequently decreasing 
by a factor of ~ 1 .5 to z ~ 0.1. 

4. SUMMARY AND CONCLUSIONS 

In this paper, we have measured the rest-frame V-band LFs 
of galaxies at 0.4 < z < 4.0 from a composite NIR-selected 
sample constructed with the deep NMBS, the FIREWORKS, 
the very deep FIRES, and the ultra-deep HST optical and 
NIR data over the HUDF, HDFN, and GOODS-CDFS fields, 
all having high-quality optical-to-MIR data. This sample is 
unique as it combines, in a self-consistent way, the advantages 
of deep, pencil beam surveys with those of shallow, wider 
surveys, allowing us to (1) minimize the uncertainties due to 
cosmic variance, and empirically quantify its contribution to 
the total error budget; and (2) simultaneously constrain the 
bright and faint ends of the rest-frame optical LFs with un- 
precedented accuracy over the entire targeted redshift range. 

The main results of our derivation of the rest-frame V-band 
LF of galaxies from z = 4.0 are (1) the faint-end slope remains 
fairly flat and does not seem to evolve from z = 4 to z = 0, 
consistent with a = -1.27 ±0.05; (2) the characteristic mag- 
nitude has dimmed by 1 mag from M v = -22.8 at z ~ 3.7 to 
My = -21.8 atz~ 0.6, and by - 1.3 mag to z = 0.1; (3) the 
normalization $* has increased by a factor of ~ 6 fromz ~ 3.7 
to z ~ 0.6, and by a factor of ~ 8 to z = 0. 1 ; and (4) the lumi- 
nosity density peaks at z ~ 1-1.5, increasing by a factor of 
~ 4 from z = 4.0 to z ~ 1-1.5, and subsequently decreasing 
by a factor of ~ 1 .5 by z = 0. 1 . 

In contrast with the findings of iRvan et al.l d2007l) . we find 
weak, if any, evidence for a steepening of the faint-end slope 
with redshift out to z = 4 . The top panel o f Figure|5]shows the 
best-fit model derived in lRvan et alj d2007l) (dot-dashed line) . 
Using the simple linear fit a(z) = a + bz, IRvan et all d2007l) 
found b rj -0. 12, whereas we find b = -0.02 ± 0.03, consistent 
with no evolution with redshift of the rest-frame V-band faint- 
end slope over the redshift range 0.4 <z< 4.0. We note that 
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the evolution of th e faint-end slope with redshift presented in 
iRvan et all (120071) was found by combining measurements of 
the faint-end slopes in the rest-frame B and UV bands. When 
only the rest-frame B-band values are considered, the evolu- 
tion of the faint-end slope with redshift becomes marginal. 



While the rest-frame B and V bands are close enough that 
significant differences should not be expected, we cannot ex- 
clude systematic differences in the (evolution) of the faint-end 
slopes in different rest-frame optical bands. It should be fi- 
nally noted that a steepening in the faint-end slope at z > 2.5 
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might still be present, due to the large uncertainties on the 
faint-end slope at z > 3. 

The hierarchical formation scenario states that many 
dwarf galaxies at high redshifts will undergo merging over 
time. This implies a steeper faint-end slope at higher red- 
shifts. Using semi-a nalytical modeling of galaxy formation, 
iKhochfar et alJ (120071) investigated the evolution with redshift 
of the rest-frame B-band faint-end slope of the LF, finding 



lOesch et al.ll 2010at lBouwens et al.ll201 UlLee et al.ll201 ll) and 
at z < 4 (a ~ -1.6T|Arnouts et alJ 120051: iReddy & Steidell 



b ~ —Q .l, in good agreement with the results from Ry an et al 



( 20071) . The predicted evolution of a from Khochfaretal 



(2007) is also plotted in Figure [5] (purple line). Our measure- 
ments are consistent with no evolution of the faint-end slope 
over the last 12 Gyr of cosmic history, in apparent contrast 
with what is expected in the current paradigm of structure for- 
mation. Supernova feedback is thought to be mainly respon- 
sible for a flatter faint-end slope compared to the dark mat- 
ter ha lo mass slope «dm at any given redshift (iDekel & Silkl 
1986), with a DM flattening with c osmic time, as s hown by 
ACDM numerical simulations (e.g jReed et al.ll2~003l) . A con- 
stant a from z = 4 would then imply either no evolution in cudm 
in the last 12 Gyr, in contrast with ACDM simulations, or a 
redshift-dependent efficiency of stellar feedback, with super- 
nova feedback being more efficient at higher redshifts when 
«dm is also steeper. Alternatively, the non-evolving rest- 
frame V-band faint-end slope out to z = 4 could be explained 
while maintaining a constant stellar feedback efficiency, if we 
assumed that the galaxies at the faint end are hosted by dif- 
ferent dark matter halos at different redshifts. In other words, 
the observed faint galaxies at low redshifts are hosted by more 
massive dark matter halos compared to similarly faint galax- 
ies at high redshifts. Additionally, systematic effects could be 
introduced by the sampled environments of faint galaxies, if 
these environments vary systematically at d ifferent redshifts. 
In fact, as shown by Khochfa r et alJ (120071) . a is steeper and 
evolves more strongly in field environments than in cluster en- 
vironments. Consequently, if the faint-end of the LF at high 
redshifts is observationally dominated by present-day cluster 
members, the faint-end slope of the overall galaxy population 
at high redshifts would then be flatter than that for the field 
LF alone at the same redshift, potentially resulting in the ob- 
served little or no evolution with redshift of a. 

Whereas we find a fairly flat and constant faint-end slope 
from z = 4 in the rest-frame V band, significantly steeper 
slopes are consistentl y found in the rest-frame far-UV, both at 
Z > 4 (a < -1.7; e.g jBouwens et al.l2007tlOuchi et al.ll2009b 



I2009t lHathi et al.ll2010t lOesch et al.ll2010bl) . Our NIR com- 
posite sample will allow us to study the evolution of the LFs 
from z = 4 at different rest-frame wavelengths, from the far- 
UV to the NIR, and to investigate the wavelength-dependence 
of the faint-end slope and its potentially different varying evo- 
lution with redshift in different rest-frame wavebands. 

We stress that we have measured the observed rest-frame V- 
band LFs over the targeted redshift range 0.4 < z < 4.0. No 
attempt was made to infer the intrinsic, dust-corrected, rest- 
frame V-band LF of galaxies, which is beyond the scope of 
our work. We note however th at, similar to what has been 
found in the local universe (e.g., iDriver et al.ll2007h . internal 
dust extinction is likely to play a significant role in shaping 
the observed LF, especially in the early universe where larger 
amount of dust extinction is typically found, even in galaxies 
around the character istic mass (e.g., iMarchesini et al.ir2010t 
lBrammeret"aIll2011l) . 

The ongo i ng Hubble Space Telescop e CANDELS 
(Grogi n et all 1201 It iKoekemoer et all 1201 ll) and 3D-HST 
(van Dokkum et al.l201 it Brammer et al.l2012l) surveys, once 
completed, will result in significantly improved constraints at 
the faint end, providing 1) very deep data over an area larger 
by a factor > 10 with respect to the area surveyed by the ultra 
deep HST data used in this work, and 2) spectroscopic red- 
shifts for thousands of faint galaxies, significantly reducing 
the uncertainties due to photometric redshift errors. 
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